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We investigate the phase structure of three-flavor QCD in the presence of finite quark 
chemical potential p/T < 1.2 by using the non-perturbatively 0{a) improved Wilson fermion 
action on lattices with a fixed temporal extent Nt = 6 and varied spatial linear extents 
Ns = 8,10,12. Especially, we focus on locating the critical end point that characterizes the 
phase structure, and extracting the curvature of the critical line on the plane. For 

Wilson-type fermions, the correspondence between bare parameters and physical parameters 
is indirect. Hence we present a strategy to transfer the bare parameter phase structure to 
the physical one, in order to obtain the curvature. Our conclusion is that the curvature is 
positive. This implies that, if one starts from a quark mass in the region of crossover at zero 
chemical potential, one would encounter a first-order phase transition when one raises the 
chemical potential. 
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I. INTRODUCTION 

At zero baryon number density, on the two-dimensional plane spanned by the light (up-down 
degenerated) quark mass mud and strange quark mass the first order phase transition around 
the massless point mud = = 0 becomes weaker as the quark masses increase, and eventually 

turns into a crossover at some finite quark masses. The boundary between the first order phase 
transition region and the crossover region forms a line of second order phase transition, called the 
critical end line. 

A question of obvious importance is the location of the critical line. Monte Carlo results on 
this issue are rather confusing at present. For the staggered fermion action, recent studies with 
improved action could place only an upper bound on the three-flavor degenerate critical quark mass, 
m, which is very small in the range of ~ 0.1 . This is in contrast to recent as well as 

an earlier study with the naive action [3-[^ which observed first order signals up to ~ 2—3. 

Furthermore, our recent study with the Wilson-clover fermion action [7|, motivated in part by the 
unclear status with the staggered action, could identify the critical end point, although the cut-off 
dependence of the location is rather large. 

The location of the critical end point in the QCD phase diagram with finite density is also an 
important issue. The first serious study on critical end point in QCD was given by Fodor and 
Katz who employ Lee-Yang zero analysis . After this study, various attempts were made and 
here we quote some reviews EQ about such studies. In this article, we address an issue of 
how the critical end line extends when switching on the chemical potential. An interesting result 
was reported in BQ which explored the imaginary chemical potential approach with the naive 
staggered fermion action. There it was observed that the critical surface has a negative curvature in 
the /i direction. This means that a first-order phase transition at zero chemical potential disappears 
when the chemical potential is increased, rather contrary to one’s naive guess. Our purpose in this 
paper is to study this question by simulations with real chemical potential using the Wilson-clover 
fermion action. This is a natural sequel of our work in B. 

The rest of the paper is organized as follows. In section m we explain a strategy on how to 
draw the critical line on the fjL-m-,^ plane. Simulation details including the parameters and the 
simulation algorithm are summarized in section IIIIl We present numerical results in section nYi 
Finally, concluding remarks are given in section IVl 


II. STRATEGY 

Let us explain our strategy to survey the phase space for Af = 3 QCD in order to identify the 
critical end point for the Wilson-type fermions. The final goal of this section is to show how to 
obtain the curvature of the critical end line on the n-m-^ plane. Note that in this section we do 
not use lattice units when expressing dimensionful physical quantities. 

First we consider the zero density case. Since the quark masses are all degenerate, we have only 
two bare parameters j3 and k {ajj. = 0 plane in the left panel in Figured]). For a given temporal 
lattice size, say At = 4, by using the peak position of susceptibility or zero of skewness of quark 
condensate, one can draw the line of finite temperature transition (the solid red line and the dotted 
green line in the left panel in Figure [T|)- The transition changes from being of first order to cross 
over at a second order critical end point (the blue point in the left panel in Figure [T|)- We compute 
the kurtosis (which is the Binder cumulant minus three) of quark condensate along the transition 
line for a set of spatial volumes A^. The intersection point is identified as the critical end point 
B. In this way, we can determine the critical end point in the bare parameter space (/3e, ke) and 
this procedure can be repeated for other values of At. 
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FIG. 1. Strategy: The left panel is the phase diagram for bare parameters spanned by j3, k and a/i for 
iVf = 3. The right panel is the same phase diagram but depicted for physical parameters spanned by mps/A 
and /i/A where A is some reference physical quantity at zero density. The blue line extending from the 
critical end point at a/i = 0 is the critical line. We study the signature of the curvature of the critical line 
with fixed W = 6. 


In order to translate the critical end point in the bare parameter space to that in the physical 
parameter space, we measure dimensionless ratios of pseudo-scalar meson mass and some reference 
quantity with mass-dimension one mps/A for the bare parameters (/3e, ke) by a zero temperature 
simulation. One can choose any reference quantity A, say T (temperature), 1/y/tQ (Wilson flow) 
0 or my (vector meson mass). To avoid the multiplicative renormalization issue, we use mps in 
the numerator of the ratio and not quark masses. In this way we pin down the critical end point 
(the blue point in the right panel in Figure [1]) in the physical parameter space whose axes are given 
by mps/A and /r/A. By repeating the same calculation for increasingly larger values of W, we can 
take the continuum limit (the orange downward arrow in the right panel in Figured]) of the critical 
end point in the physical parameter space at zero density. 


^p°S^e(/^ — 0) _ ,. "rps,E(l^ = 0) 

A|?“t(/r = 0) “ Nt^oc Ae(/U = 0) ^ ’ 

This strategy is in fact used in our zero density study [^. 

When switching on the chemical potential, the basic procedure is the same; one just has to 
repeat the same analysis on a different plane with /r 7 ^ 0 (See the left panel in Figure [T])- For a 
fixed lattice temporal size, Nt = 6 , in order to draw the critical end line, we consider a pair of 
dimensionless ratios 


"lPS,E(/.i) , /.i 

mps,E(0) '''' rE(o)’ 


( 2 ) 


where for each ratio we have chosen proper reference quantities at zero density. By plotting 
these two quantities one can obtain a critical line as shown in the right panel of Figure dl We 
are interested in seeing whether the critical line bends toward the lighter mass or heavier mass 
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direction. More quantitatively, from a fitting 

"1ps,e( 0)J ^ V7rrE(0)J ^V7rrE(0)J 


(3) 


we shall extract the curvature ai and see its sign, and this is the final goal of this paper. 

If one wants to take the continuum limit of the critical end line, one has to take the Nt ^ oo 
limit for fixed values of fi/T^iO) 


mP>s°E(M) 


"iPS,E(M) 


lim 

Nt^oo mpg^E(O) 


fixed At/Tsfo) 


(4) 


After repeating the same procedure with different values of fi/T^iO), one can plot (?t^p°^e(i^)/^p°s^M 0))^ 
as a function of /u/Te(0). Then by fitting with the same form as in eq.([3]), one can obtain the 
curvature in the continuum limit. 


III. SIMULATION DETAILS 


We employ the Wilson-clover fermion action with non-perturbatively tuned Cgw iM ™ 
presence of chemical potential with the anti-periodic boundary condition in the temporal direction 
for fermion fields while the periodic boundary condition is imposed for spatial direction. The 
Iwasaki gauge action [1^ is used for the gluon sector and gauge link variables satisfy the periodic 
boundary condition. The number of flavor is three, Af = 3, and the masses and chemical potentials 
for quarks are all degenerate. The temporal lattice size and the simulated quark chemical potential 
are fixed to At = 6 and a// = 0.1, respectively, and thus /i/T = 0.6. In our study, the phase 
reweighting method explained below is used to deal with the complex phase, and to survey a wide 
range of /x and k, we adopt the multi-parameter reweighting method; details are given in Appendix 
El To perform finite size scaling analysis, the spatial volume is changed over the linear sizes Ag = 8, 
10 and 12. In order to search for the transition point, we select four /3 points (/3 = 1.70, 1.73, 1.75 
and 1.77) and for each /3, we vary k to locate the transition point. 

The phase reweighting method is adopted to handle the complex phase according to 


{ 0 ) = 




(5) 


where (•■•)|| is the average with phase quenched fermion determinant 




( 6 ) 


and the phase factor for one-flavor is given by 




det D(iu) 

I det D{^)\ 


(7) 


Configurations are generated by RHMC [16| | with the phase quenched quark determinant. The 
MD step size is chosen such that a reasonable acceptance rate > 80% is retained. For each lattice 
parameter set (/3 ,k, At, Ag) we generate 0(100,000) trajectories and the configurations are stored 
at every 10th trajectory; the order of number of configurations are 0(10,000) for each parameter 
set. The phase factor and //-derivatives of the fermion determinants required in //-parameter 
reweighting are computed exactly using the analytical reduction technique [iM3 for all stored 
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FIG. 2. The average of phase reweighting factor with iVf = 3 as a function of k. /r/T = afi x Nt = 0.1 x 6 = 
0.6. The reweighting factor is significantly away from zero. This shows that the sign problem is mild in this 
region. 


configurations. The dense matrix obtained by the reduction is numerically computed on GPGPU 
with LAPAGK routines. We measure the trace of quark propagator and its higher power up to 
fourth order which are used not only for the computation of higher moments of quark condensate 
but also for the parameter reweighting (See Appendix!^ for details). In the computation of traces, 
we adopt the noise method with 20 Gaussian noises that is checked to be sufficient to control the 
noise error. 

For each fixed parameter set (/3, a^, Nt, Ng), we make runs at several values of k. In order 
to integrate those runs we adopt the multi-ensemble reweighting technique and search for 
the transition point in k for the fixed parameter set. See Appendix |B] for the details of the multi¬ 
ensemble reweighting. Here we only mention that we use some approximation to efficiently evaluate 
the quark determinant in the reweighting factor as well as observables at many reweighting points. 

In our approach, there are practically two important issues: the overlap problem and the validity 
of approximation made at calculating the ratio of quark determinant in the reweighting factor. The 
issue of the overlap problem will be addressed in the next section. The validity of the approximation 
is discussed in Appendix |A] and the conclusion is that the approximation we made is safe in our 
parameter region. 

The physical scale settings we use in this paper, for example the Wilson flow scale ^/to and 
the hadron mass, are taken from Ref. [^. 


IV. RESULTS 
A. Phase reweighting factor 

Figure [2] shows the average value of the phase-reweighting factor as a function of k. For small k 
and large volumes, the value becomes smaller, signaling that the sign problem is becoming serious. 
Nevertheless, it stays away from zero (> 0.5) beyond statistical error, guaranteeing the validity of 
the phase-reweighting for our range of lattice parameters. 

We also check the average of the reweighting factor wt in eq. (IBOp of the multi-ensemble reweight¬ 
ing, (r(;T)ME = wt:{U)/ Yhs where unexplained notation is given in Appendix iBl The relative 
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FIG. 3. Relative error of the reweighting factor of multi-ensemble reweighting in ea. (IB3|) as a function of n 
for P = 1.70 (left-top), 1.73(right-top), 1.75(bottom-left), 1.77(bottom-right). The shaded bands represent 
the respective k values at transition/cross-over, from Table U and the dotted lines indicate the location of 
the simulated k values. Here only three selected values of the chemical potential are shown for each (3. The 
spatial lattice is fixed Ns = 12. The relative weight is significantly smaller than one, thus the sign problem 
is mild even for the reweighted parameter space. 


error of the reweighting factor is plotted in Figure [3l The errors are estimated by the jackknife 
method with bin size of 1000 configurations. Figure [3] shows that the relative error is sufficiently 
small [error of (u)T)ME]/(tfT)ME I, even at larger chemical potential afi « 0.2. This means that 
the central value of the reweighting factor is significantly away from zero beyond many sigmas. 
Thus, we conclude that the overlap problem is not so severe in our parameter region. 


B. Moments of chiral condensate and transition point 

At finite quark mass, the quark bilinear operator 'ijj'ijj is not a real order parameter but considered 
to be a mixture of “energy” and “magnetization” operators 0]. We study the bilinear operator as 
a primarily magnetization operator, however, since we do not have enough data set to resolve the 
mixing of observables. The detailed practical definition of its moments is given in Ref. j^. 

Figure 0] shows curves of the susceptibility and kurtosis for quark condensate obtained by the 
multi-ensemble reweighting. The error bands are estimated by the jackknife method with bin size 
of 500 — 1000 configurations. For = 0.1, the averages at each point of data generation are 
shown in order to illustrate how multi-ensemble curves interpolate those raw data. At /? = 1.73, 
the curves reweighted to a/r = 0 can be compared with data generated at zero density [^. The 
agreement supports the validity of multi-ensemble reweighting and jackknife error estimation away 
from a/i = 0.1. The applicable range of /i/K-reweighting depends on /3, and judged from the growth 










susceptibility 
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FIG. 4. The susceptibility and kurtosis of quark condensate as a function of k for /3 = 1.70 (left-top), 
1.73(right-top), 1.75(bottom-left), 1.77(bottom-right). In each figure, in addition to the raw data at afi = 0.1 
with Ns = 8,10,12, we plot the reweighting results expressed by a band with selected 3 values of a/r. All 
reweighting results are produced from configurations at a/i = 0.1 and several ensembles with different k 
values are integrated into the multi-ensemble reweighting. In /3 = 1.73 with Ng = 10,12, we also plot the 
raw data at a/i = 0 given in the zero density study Q- These raw data and the reweighting results are 
consistent with each other, although they are completely independent. This shows that the reweighting 
together with the approximation used for the reweighting factors and observables (See Appendix for 
details of the approximation) is fine. 


of error, the lower /3 tends to have a larger applicable range. 

As seen in the figures, the locations of the maximum of susceptibility and minimum of kurtosis 
are consistent with each other. Furthermore the skewness zero location is also consistent with 
them although it is not shown here. We take the location of the maximum of susceptibility as 
the transition point. The numerical values are summarized in Table H] where the peak height of 
susceptibility Xmax and the minimum of kurtosis iFmin are also listed for selected values of a/r. 

As seen in Table HI the volume dependence of the transition points is rather mild. Hence the 
thermodynamic limit can be safely taken with a fitting ansatz, 

utiNs) = Kt{oo)+ c/N^. 


(8) 
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FIG. 5. The phase diagram of A^f = 3 QCD with finite chemical potential projected on the (/3, k) plane. 
The transition points are expressed by open symbols while the critical end points are given by filled ones. 
The lower /3 side of the critical end point is the first order phase transition region. For larger chemical 
potential, the critical point moves toward the upper-left corner. The Kc line where the pion mass vanishes 
is also shown. 


The resulting value of Kt{oo) is shown in Table HI The phase diagram of bare parameters f3 and k 
is given in Figured The transition lines have a sensitivity on the value of chemical potential, afi. 


C. Kurtosis intersection 

The next step is to determine the critical end point. For that purpose we adopt the kurtosis 
intersection method 0]. The value of kurtosis can be used to diagnose the strength of phase 
transitions. For a first order phase transition, the infinite volume value of kurtosis is K = —2 while 
for a crossover it is iF = 0. At the critical point as the end point of a first order phase transition 
line, the kurtosis is expected to take the same value irrespective of the spatial volume between —2 
and 0. The value at the critical end point depends on the universality class of the second order 
phase transition. 

Figure [6] plots the minimum of kurtosis as a function of /3 for some selected values of a^. This 
shows that a strong first order phase transition at lower /3 becomes weaker for higher (3 and such 
a change becomes rapid for larger volumes. We fit the data with the fitting form inspired by 
finite size scaling, 


K^in = Ke + (/3 - /3e) , (9) 

where Ke, A, v and /3 e are fitting parameters and the results are listed in Table HIl The resulting 
exponent v and the value of kurtosis at the critical end point Ke are independent of afj, within 
errors, and they are consistent with the values of 3-dimensional Z 2 universality class, u = 0.63 and 
Ke = —1.396 respectively. On the other hand, the universality class of 3-dimensional 0(2) and 
3-dimensional 0(4) are rejected, rather strongly by the value of Ke- 
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FIG. 6. Kurtosis intersection at a/r = 0.00 — 0.19. In the fitting, three lowest values of f3 are used. The 
black pentagon represents the critical end point (CEP) in bare parameter space /3 e and it moves to the lower 
side for larger chemical potential. See Table [TTl for the values of fitting parameters. The horizontal magenta 
line shows = —1.396 for 3-dimensional Z 2 universality class. In this region of the chemical potential, the 
value of is constant, namely the universality class does not change. 


We superimpose the obtained critical end points (/3E(a/i), for 0 < < 0.19 in the 

phase diagram of Figure [5j The critical end point moves toward the upper-left corner by increasing 
a/i. 

In order to confirm the universality class and the location of the critical end point, we check 
another exponent y/iz which is obtained from the volume scaling of the susceptibility peak height 
of quark condensate 

Xma. = CNl 


(10) 
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j3 an Ns 

Xmax -^min 

1.70 0.00 8 

1.70 0.00 10 

1.70 0.00 12 

1.70 0.00 00 

0.1415100(52) 63.25(57) -1.6237(72) 

0.1415209(28) 121.74(91) -1.7447(50) 

0.1415177(15) 210.09(90) -1.8336(24) 

0.1415203(27) 

1.70 0.10 8 

1.70 0.10 10 

1.70 0.10 12 

1.70 0.10 00 

0.1414498(52) 60.54(57) -1.5946(83) 

0.1414577(28) 112.71(93) -1.6951(63) 

0.1414541(15) 190.9(1.0) -1.7826(38) 

0.1414554(31) 

1.70 0.20 8 

1.70 0.20 10 

1.70 0.20 12 

1.70 0.20 00 

0.1412465(59) 48.22(65) -1.424(15) 

0.1412493(32) 76.0(1.7) -1.417(25) 

0.1412516(27) 125.0(7.3) -1.550(64) 

0.1412537(46) 

1.73 0.00 8 

1.73 0.00 10 

1.73 0.00 12 

1.73 0.00 oo 

0.1404267(44) 25.66(20) -1.2943(84) 

0.1404304(37) 40.90(51) -1.310(12) 

0.1404340(16) 54.54(48) -1.2600(94) 

0.1404371(29) 

1.73 0.10 8 

1.73 0.10 10 

1.73 0.10 12 

1.73 0.10 00 

0.1403406(44) 24.68(20) -1.2580(87) 

0.1403400(38) 37.29(47) -1.228(13) 

0.1403421(16) 47.41(44) -1.126(11) 

0.1403427(30) 

1.73 0.19 8 

1.73 0.19 10 

1.73 0.19 12 

1.73 0.19 00 

0.1401016(50) 20.60(22) -1.082(19) 

0.1400964(59) 26.22(82) -0.983(53) 

0.1400946(92) 29.5(1.8) -0.72(13) 

0.1400913(97) 

1.75 0.00 8 

1.75 0.00 10 

1.75 0.00 12 

1.75 0.00 00 

0.1396591(87) 14.37(21) -1.115(16) 

0.1396684(58) 19.41(41) -0.985(21) 

0.1396682(37) 23.27(42) -0.878(20) 

0.1396722(64) 

1.75 0.10 8 

1.75 0.10 10 

1.75 0.10 12 

1.75 0.10 00 

0.1395533(87) 13.82(21) -1.077(17) 

0.1395547(59) 17.90(39) -0.914(21) 

0.1395546(39) 20.52(38) -0.765(21) 

0.1395552(67) 

1.75 0.18 8 

1.75 0.18 10 

1.75 0.18 12 

1.75 0.18 00 

0.1393077(97) 11.86(22) -0.938(28) 

0.1392962(75) 14.37(38) -0.776(42) 

0.1393035(97) 14.93(53) -0.59(10) 

0.139295(12) 

1.77 0.00 8 

1.77 0.00 10 

1.77 0.00 12 

1.77 0.00 00 

0.1389157(85) 8.34(15) -1.009(16) 

0.1389300(88) 10.17(29) -0.794(31) 

0.1389349(93) 11.41(43) -0.619(50) 

0.138944(11) 

1.77 0.10 8 

1.77 0.10 10 

1.77 0.10 12 

1.77 0.10 00 

0.1387866(88) 8.02(14) -0.976(16) 

0.1387902(98) 9.36(29) -0.708(34) 

0.138792(10) 10.34(40) -0.549(51) 

0.138794(13) 

1.77 0.16 8 

1.77 0.16 10 

1.77 0.16 12 

1.77 0.16 oo 

0.1385825(95) 7.29(14) -0.914(20) 

0.138570(13) 8.07(27) -0.610(45) 

0.138559(15) 8.98(41) -0.453(88) 

0.138553(18) 


TABLE 1. The transition point Kt, the peak hight of susceptibility and the minimum of kurtosis. The errors 
are estimated by the jackknife analysis except for the value of Kt aX, Ns = oo where the error is calculated 
from the fit in eg.dSl). 
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a/j, 


/3e 


0.00 

0.05 

0.10 

0.15 

0.16 

0.17 

0.18 

0.19 


1.7247(10) 

1.72249(92) 

1.71803(76) 

1.71372(95) 

1.7129(11) 

1.7119(14) 

1.7109(20) 

1.7095(32) 


ke_ 

0.1406320(29) 

0.1406904(29) 

0.1407928(29) 

0.1408541(31) 

0.1408622(33) 

0.1408723(38) 

0.1408864(48) 

0.1409120(62) 


Ke _ 

-1.366(15) 

-1.383(14) 

-1.403(12) 

-1.381(15) 

-1.371(18) 

-1.358(22) 

-1.346(28) 

-1.335(43) 


0.734(35) 

0.694(32) 

0.615(27) 

0.569(33) 

0.567(38) 

0.570(49) 

0.590(74) 

0.66(13) 


A _ 

0.64(10) 

0.545(86) 

0.371(64) 

0.281(67) 

0.279(78) 

0.285(99) 

0.33(15) 

0.47(33) 


X^/dof 

3.76 

3.24 

3.18 

5.54 

4.94 

3.89 

2.86 

2.03 


Universality class 

Ke 

V 

itv 

3-dimensional Z 2 

-1.396 

0.63 

1.964 

3-dimensional 0(2) 

-1.758 

0.672 

1.962 

3-dimensional 0(4) 

-1.908 

0.748 

1.975 


TABLE 11. Fit results for kurtosis intersection for selected values of a/i. The errors are estimated by the 
jackknife method, x^/dof is the average value. In the region 0 < a/i < 0.19, the exponent and the value of 
kurtosis at the critical end point Ke are constant within errors and the universality class is consistent with 
3-dimensional Z 2 , while other universality classes, 3-dimensional 0(2) and 3-dimensional 0(4) are rejected. 


with fit parameters b and C. The exponent b depends on the nature of transition, i.e., b = d the 
spatial dimensionality at a first order phase transition, and 6 = 0 for a crossover. At the critical 
point as the boundary of the first order phase transition line, the exponent is expected to be 
6 = 7 / 1 / with critical exponents 7 and z/. Figure [7] shows the exponent 6 along the transition line as 
a function of /?. We observe that the exponent 6 at the critical end point estimated by the kurtosis 
intersection is consistent with the value for the 3-dimensional Z 2 , 7 / 1 / = 1.964. Thus we observe a 
consistency between the kurtosis intersection analysis and the volume scaling of susceptibility. We 
note that it is difficult to differentiate universality classes depending solely on the volume scaling 
of the susceptibility peak since the values of 7 / 1 / are very close to each other as listed in Table HH 


D. Critical line 


The analysis of the critical line below requires a careful manipulation with scale setting. Thus 
we distinguish quantities in lattice units from those in physical units by placing a tilde on the 
former, e.g., the chemical potential in physical units is denoted as /i and that in lattice units by 
/i = a/i. 

In the previous subsection, we have determined the critical end points in the bare parameter 
space. The last step is to translate the critical end point on the (/3, k) plane to the physical space, 
to obtain the critical line as one varies /i, and finally to extract its curvature. For that purpose, as 
explained in Sect. m we need to compute the pair of ratios in eq.([ 2 ]) as follows. 


wips,E(7t) _ dT-ps,E(/t) a(0) 
”^ps,e( 0) mps,E(0) a(/i)’ 
/i _ a(0) 

tWoT ^ ^ 
rE(0) a(/i) 


( 11 ) 

( 12 ) 


where mps,E(/l) is the pseudo-scalar (PS) meson mass in lattice units evaluated at (/3, k) = 
(/3e(/1))'^e(/ 1))- Note that the PS mass at the critical point does not depend on /t directly, but 
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FIG. 7. Exponent b of peak height of susceptibility as a function of /3 (open symbols). The lines connecting 
points are just for a guide for your eyes. The filled regions represent the corresponding critical end point 
determined by the kurtosis intersection method. The horizontal three lines {b ~ 2) represent the values of 
7 / 1 / for universality classes Z 2 , 0(2) and 0(4) in 3-dimension. It is hard to see their difference at this scale. 


only through /3 e and ke at /2. The PS mass is measured by the zero temperature simulation at /3 e 
and ke- On the other hand, the lattice spacing requires some careful thought as follows. 

We usually determine the lattice spacing by choosing a line of constant physics (TCP) and 
specifying the value of a dimensionful physical quantity on that line. For example, one may choose 
the dimensionless combination mps\/to for specifying the TCP, and the value of mps in physical 
units to determine the lattice spacing along the chosen TCP, 


a{l3,y) 


mpsif3,Ky{/3)) 

mps{y) 


(13) 


where y is the value of the constant physics y = mps^/to and Ky(/3) is defined such that the 
following equation holds for each /3, 


y = hrps(/3, Ky{/3))\^{(3, Ky{l3)). (14) 

The notation of the lattice spacing in ea. (IIIII2p means that 

«(A) = «(/^E(h),y)- (15) 

Note that, again, the lattice spacing does not depend on y directly, but only though the /3e at y. 
Thanks to TCP, where the physical unit mass in the denominator in eq. (|I3p is not known a priori 
but common, the physical mass cancels out in the ratio of lattice spacings and the ratio may be 
computed by using the PS mass in lattice units, 

a(0) ^ mps(/3E(0),Ky(/3E(0))) 

a(y) 'rhps{l3Eiy),HyiPE{y)))' 
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constant physics 

scale input 

ai 

02 

uips-Zh) = 0.55 

1 / 

1.924(60) 

-0.58(72) 

mps\/to = 0.65 

1 /-/^ 

2.148(39) 

-1.74(52) 


TABLE III. The curvature of the critical line in the fitting form in eq.([3]) for scale inputs (i/to) and the 
values of constant physics TOps\/to = 0.55,0.65. The error of curvature is estimated by the jackknife method. 


In the following, for the computation of the ratio of lattice spacings, we use the Wilson flow scale 
instead of the PS mass since the former is precisely calculated 

«(0) ^ 1/V^(/e(0), k^(/3e(0))) 

One can employ a different LCP by specifying a different value of y// y). The resulting 
lattice spacing coincides with that from the original (y) definition if, in specifying the value of the 
dimensionful quantity, one takes into account the variation of that quantity in moving from the 
original LCP to a new LCP. In general the agreement will not be exact due to scaling violations. 

a{l3,y') = a{j3,y) + (lattice artifacts). (18) 

Thus differences one may observe in physical results due to the choice of LCP is a scaling violation 
effect. In the following, we choose two values for the line of constant physics, 

y = mpsVto = 0.55 and 0.65. (19) 

We use the Wilson flow scale and the hadron mass computed in Appendix A of Ref. , where 
the zero temperature simulations were carried out with the same lattice actions and sufficiently 
large lattices > 5. Especially, we select /3 = 1.70, 1.73, 1.75 and 1.77 data in our analysis 
here. By combining the above scale inputs and the information of the critical end point at finite 
chemical potential determined in the previous subsection, we calculate the two ratios in eq. m 
and (fT^ . The results are plotted in Figure [8l 

We extract the curvature by using the fitting form in eq. ([3]) . The results are tabulated in Table 
imi The errors of fitted parameters are estimated by the jackknife method using the uncorrelated 
chi squared function in each ht. We also try to perform a fit including correlations by using the 
covariance matrix estimated by the jackknife method; the results are consistent with the above 
analysis although the covariance matrix is poorly estimated. We observe that the critical line 
has a sensitivity on the value of constant physics. This difference is considered as a systematic 
uncertainty caused by the choice of the scale setting as discussed above. All in all, we hnd the 
curvature of the critical line to be positive with a statistical error of about 3% and a systematic 
error of about 10%. 


V. CONCLUDING REMARKS 

We have investigated the critical line on the y-ruTr plane, especially its curvature, in Af = 3 
QCD by using non-perturbatively 0{a) improved Wilson fermion action. We have determined the 
critical end point by making use of the kurtosis intersection method. The critical line is drawn by 
repeating the calculation in the range of chemical potential with applications of various reweighting 
techniques, that is, the multi-parameter/phase/multi-ensemble reweighting. 
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^/{7iTe(0)) 


1 /2 

FIG. 8. Critical line for constant physics mps^o^ = 0.55,0.65 with scale input: the Wilson flow scale. 
Both lines are bending toward the upper direction. 


The value of kurtosis at the critical end point and the exponent v obtained from the kurtosis 
intersection analysis in the range of chemical potential we investigated are consistent with those of 
the 3-dimensional Z 2 universality class. Furthermore, if the above universality class is used as an 
input in the analysis of exponent extracted from susceptibility peak, the expected location of the 
critical point is consistent with that obtained from the kurtosis intersection method. 

Our analysis shows that the curvature of the critical line is positive. This disagrees with a 
previous study with the naive staggered fermion action [^, 12] where the critical line is expressed in 


terms of quark mass. We note that neither the previous study nor ours have taken the continuum 
limit. Thus further work with larger W is desired. 
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Appendix A: Reweighting 


The phase reweighting and multi-parameter reweighting for bare parameter^ mo and ^ can be 


done by the formula 


(C>(mo,/r')) 




del D{mo,fi) J 






det J 




(Al) 


Here, ttiq and /r' are target parameters while mo and /r are actual simulation parameters. The 
average ™ LHS is taken by using the Boltzmann factor including the full quark determinant 

at parameter mg, /r', while the average (•••)||,mo,At ™ RHS is taken by using the Boltzmann factor 
including the phase quenched quark determinant at parameter mo, R- Here we have explicitly 
written down the bare parameter dependence on the observable O(mo, ^), say the quark propagator. 

In ea. (|Al|l . one needs to evaluate the reweighting factor. 


det D{mo, fj,) J ’ 


(A2) 


where the second factor is already computed and stored but the first factor, the ratio of quark 
determinants, requires high cost computation if one tries to calculate it directly at many target 
parameter points (mg,//'). Thus we adopt a cheaper approximation method, that is, the Taylor 
expansion of the logarithm of determinant which is known to have better convergence property 
than the other expansion schemes (^ . 


In 


/det D{mQ, /r') 
\ det D{mo, iJ.) 



j'-kl 




k 

In det D{mo, /r) 


—Indet T)(mo,/i), (A3) 


^ We use the bare mass parameter mo instead of ft = 1/ (2mo-|-8), since the former parameter is useful in the following 
discussion. In this appendix, fj, is the chemical potential in lattice units. We do not consider /3-reweighting. 
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Amo = iTT'O ~ ^0) (A4) 

— fi. (A5) 


Once some leading coefficients in the expansion are calculated, one can easily evaluate the ratio 
at many reweighted points up to truncation errors. In our calculation, we include the following 
coefficients 


(j> k) 


(1 — 4,0) : purely mo-derivatives, 

(0,1 — 4) : purely //-derivatives. 


(A6) 


but no mixed derivatives say (1,1) and so on. The explicit form of the approximated ratio of 
determinant is given by 


/det D{mQ, //') 
\ det D{mo, /t) 


TVf 


exp 


Ni //) + N^Y Wk{mo, fi/T) 

j=i ^ k=i 


(A7) 


where Vbi, 2 , 3,4 are given in [22l |. 

For the same reason as the ratio of determinant, we use an expansion form of the observable in 
eq. flAT]) 


O{m'o,n') = Y 

j,k=0 


Am;^A/i^ 

f-kl 




k 

O{mo,^). 


(A8) 


For the trace of higher powers of quark propagator which is included in the higher moments of 
quark condensat^, we apply the following approximation (//-derivative terms are totally neglected) 


3 


trT)“^(mo,//') ? 

atrD ^(mo,//)-k ^(—)-^ Am;jtri7 
i=i 

(A9) 

trF>“^(mo,//') ? 

2 

a trF>"2(mo,//) -k ^(-)-^(j -k l)Am^tri:»"(-^+2)(mo,//), 
j=i 

(AlO) 

trT)“3(m'o,//') ? 

1 

a trF>“3(mo,//) -k ^(-)-^(j -k l)(j -k 2)Am^trT)"(^+^)(mo,//), 

7 = 1 

(All) 

trT)“^(mo,//') ? 

J 

trF)“'^(mo,//). 

(A12) 


At first glance, you may doubt the approximation for the ratio of determinant in ea. (IA7h and 
observables in ea. (|A9tiA12p especially the higher powers of the inverse. We however have some 
evidences that this approximation is good within our parameter range and statistical error. First, 
the approximation for the observables in eq. ()A9IIA12]) is compared with the partial quenching results 
where there is no truncation error in the observable. Even for the kurtosis including ea. (IA12P . we 
do not see any significant difference between them within errors. Second, in order to check the 
effects of the mixed derivative terms in the reweighting factor, we compute and include the mixed- 
derivative coefficients up to 4th order, namely {j,k) = (1,1), (2,1), (3,1), (1,2), (2,2), (1,3) 
coefficients in ea. (IA3p . and check their effects on the moments of chiral condensate (of course. 


^ The formulae for higher moments of quark condensate in terms of quark propagator are explicitly given in Ref. [3] . 
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the associated mixed derivative contributions for the observable are also included) at /3 = 1.73 
and Ng = 12 in the range of the chemical potential, 0 < /r < 0.19. Then it turns out that the 
difference is quite small, that is, within statistical errors in the parameter space. Furthermore, we 
check the hierarchy of the terms and observe that the dominant term is (1,0) and the leading term 
in the mixed derivative terms is (1,1), and then it turns out that the magnitude of their ratio, 
|[(1,1) term]/[(l,0) term]| is of order 10“^ at a maximum. This shows that neglecting the mixed 
derivative terms is justified and our approximation in eq. (IA7h is fine. We naturally expect that the 
same goes for the other cases of /3 = 1.70,1.75,1.77. Thus, we conclude that the approximation 
made in eq. (IA7p is legitimate in the range of parameter we explored. 


Appendix B: Multi-ensemble reweighting 


In this appendix, we review the multi-ensemble reweighting technique in 20|. An estimated (E) 
expectation value of some operator at a target parameter denoted by T is given by 


oE EuWT{U)nT{U) 

^ Eu'^TiU’) ■ 


(Bl) 


where Ylu abbreviation of sum over all configurations, namely sum over all ensemble (each 
ensemble is numbered by r) and all configurations (numbered by n) therein. 


R Nr 




(B2) 


u 


r=l n=l 


with the number of ensembles R and the total number of configurations W for ensemble r. 
The reweighting factor wt is given by 




wt{U) = 


Ef=i Ns exp St{U) - S'J(U) + fs-F t 


(B3) 


where St is the action at the target parameter, and SJ is the simulated actions (using phase 
quenched determinant at the simulated parameter) for ensemble s {s = 1,2, In our cas^ 


g-ST(Wn) ^ g-SG(Wn) det D{kt, /xt; Ur,n), 

^-shUr,n) -^ D{Ks, IJ.s',Ur,n)\- 

The ratio of Boltzmann weight in eq. ()B3jl is given by 


(B4) 

(B5) 


exp 


ST{Ur,n) - SUUr,n) 


I det D(^Ks , /is, f7r,n) 1 

detD{KT,fiT; Ur,n) 

I det D{Ks-,^s\Ur,7i)\ —i9{Kr,Ur\Ur.‘n.) 

I det D{Kr,llr]Ur,n)\ 

det Z)(KT,/iT;t^r,n) 
det D ) ^-^v^TL ) 


(B6) 


where we have already measured the phase 9{Kr, fJ,r,Ur,n)- The ratio of determinants can be 
estimated by using the expansion method given in the previous appendix. 


® In contrast to the previous appendix, we use k instead of the bare mass mo in this appendix. The arguments of the 
Dirac operator are k, and gauge configuration of U. Note that r denoting the ensemble of Ur,n and s denoting 
the ensemble of parameter of action S'!' are independent. 
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fs (s = 1,2, ...,R) in ea. (IB3ll which are free parameter and we determine them by solving the 
non-linear equation. 


fs = Fs-lnY^ 


u 



exp 




(B7) 


where Fg are dummy to avoid numerical instability. We solve the equation by iteratively substi¬ 
tuting trial values of fs with initial values /^ = 0 for all s. We observe that this iteration converges 
after around (or less than) 20 iterations for all cases. 

Ft in eq. fIBOp is a constant to avoid numerical instability. 


Fp 




U s=l 


(B8) 






